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Abstract 

In  this  report  we  address  the  problem  of  nonlinear  filtering  in  the  presence  of  integer  uncer¬ 
tainty.  In  the  simulation  results  we  show  that  particle  filtering  is  capable  of  resolving  integer 
ambiguity  in  the  given  nonlinear  setup.  Motivated  by  these  results  we  introduce  particle  filtering 
for  an  exponential  family  of  densities.  We  prove  that  under  certain  conditions  the  approximated 
conditional  density  converges  to  the  true  conditional  density.  For  the  case  where  the  conditional 
density  does  not  lie  in  an  exponential  family  but  stays  close  to  it,  we  show  that  under  certain 
assumptions  the  error  of  the  estimate  given  by  this  approximate  nonlinear  filtering,  projection 
particle  filtering ,  is  bounded.  In  the  simulation  results  we  show  the  application  of  particle 
filtering  to  Global  Position  System  (GPS). 


1  Introduction 

GPS  provides  world  wide  positioning  with  acceptable  accuracy,  if  four  or  more  satellites  are  in  view. 
Although  the  satellite  constellation  guarantees  availability  of  four  or  more  satellites  (sometimes 
even  nine)  world  wide,  natural  or  man-made  obstacles  can  easily  block  the  satellites’  signal.  To 
overcome  this  vulnerability,  one  might  think  of  integrating  dead  reckoning  or  Inertial  Navigation 
System  (INS)  with  the  GPS  [1]  [2]  [3] .  In  this  case,  the  INS  or  the  dead  reckoning  provide  positioning 
that  is  adjusted  by  the  GPS. 

Using  differential  GPS  allows  the  user  to  have  a  more  accurate  measurement.  In  fact,  a  good 
portion  of  the  positioning  error  can  be  removed  from  the  estimation.  This  and  new  technology  allow 
the  use  of  the  carrier  phase  as  part  of  the  positioning  information.  This  can  increase  the  accuracy 
of  the  estimation  to  centimeter,  or  in  the  static  case,  to  millimeter  levels.  This  can  happen  only  if 
we  are  able  to  estimate  the  number  of  full  cycles  of  the  carrier  phase,  which  cannot  be  measured. 
This  problem  is  called  integer  ambiguity  resolution  [4]  [5]  [6]  [7] . 

Although  carrier  phase  differential  GPS  allows  for  very  accurate  positioning,  it  is  very  sensitive 
to  obstacles  that  can  block  satellite  signals  and  cycle  slips.  A  good  estimation  algorithm  should  be 
able  to  quickly  estimate  the  integer  ambiguity  on  the  fly.  Most  of  the  algorithms  use  integer  least 
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square  methods  for  this  [6]  [7]  [4],  In  [4]  a  Kalman  filter  type  of  setup  is  used  to  estimate  the 
integer  ambiguity. 

In  most  of  the  applications,  integrated  INS/GPS,  dead  reckoning/ GPS,  or  vehicle  dynarn- 
ics/GPS,  linearization  of  the  dynamics  and  the  GPS  observation  is  the  main  tool  for  estimation 
[6]  [7]  [4].  It  can  be  shown  [3]  that  when  the  number  of  satellites  is  below  a  certain  level  or  the 
geometry  of  the  current  constellation  is  near  singular,  the  linearization  causes  the  system  to  be 
unobservable.  In  this  case,  it  is  important  to  use  a  nonlinear  setup  for  the  estimation  problem.  In 
[3]  this  case  is  studied  by  using  an  approximation  for  nonlinear  filtering  [8]  [9] . 

Except  for  very  special  cases  in  nonlinear  settings,  estimating  the  state  given  the  observations 
results  in  an  infinite  dimensional  filter  [10] .  Therefore,  approximation  methods  of  finite  dimension 
are  very  appealing.  The  most  widely  used  approximation  filtering  method  is  the  Extended  Kalman 
Filter  (EKF),  which  is  a  heuristic  method  based  on  the  linearization  of  the  state  dynamics  and 
observation  near  the  nominal  path  [10].  EKF  is  computationally  simple  but,  the  convergence  of 
the  conditional  distribution  to  the  actual  distribution  is  not  guaranteed. 

Projection  filtering  is  another  approximation  method  [11]  [12]  [13]  [14].  In  projection  filtering  it  is 
assumed  that  the  conditional  density  of  the  state  of  the  system  can  be  approximated  by  a  member 
of  parametric  family  of  densities.  In  this  case,  estimating  the  conditional  density  is  equivalent  to 
estimating  the  parameter  of  the  family.  In  [11]  [12]  [13]  the  exponential  family  of  densities  is  chosen 
as  the  parametric  family.  In  [14]  the  approach  is  different;  there  a  Galerkin  approximation  is  used 
for  solving  the  Fokker-Planck  equation  [10]. 

An  entirely  different  approach  to  approximate  the  conditional  density  was  proposed  in  [8]  [9]. 
This  method  is  based  on  the  Monte  Carlo  method  and  is  called  particle  filtering.  In  this  method, 
the  particles  at  time  tt  are  i.i.d.  random  vectors  that  are  distributed  according  to  the  empirical 
conditional  distribution  of  the  state,  given  the  observations  up  to  time  f.  These  particle/state 
vectors  are  used  in  the  state  equation  to  find  the  values  of  particles  at  time  fj+i.  Then  at  time 
ij+i,  the  empirical  distribution  is  evaluated  according  to  the  values  of  the  particles.  The  new 
observation  at  time  tj+i  is  taken  into  account  through  Bayes’  Rule  to  calculate  the  conditional 
empirical  distribution,  and  this  process  can  be  repeated.  In  [8]  it  is  proved  that  by  a  large  enough 
number  of  particles,  one  can  get  an  approximate  conditional  distribution  that  is  arbitrarily  close 
to  the  true  conditional  distribution. 

In  the  cases  where  we  have  some  prior  information  about  the  distribution,  we  should  expect 
to  achieve  higher  performance  if  we  take  this  information  into  account.  By  higher  performance, 
we  mean  a  reduction  in  the  computational  cost  and  an  increase  in  the  convergence  rate.  Here 
we  assume  that  the  conditional  distribution  has  a  density  in  an  exponential  family  of  densities, 
or  at  least  stays  close  to  it  in  a  sense  that  we  will  define.  Using  this  assumption,  we  replace  the 
empirical  distribution  in  [8]  with  the  Maximum  Likelihood  Estimate  (MLE)  of  the  parameters  of  an 
exponential  density.  We  call  this  new  method  projection  particle  filtering.  In  Theorem  5.6  we  show 
that  if  the  conditional  density  of  the  state  given  the  observations  lies  in  an  exponential  family  of 
densities  then  the  estimated  conditional  density  converges  to  the  true  conditional  density  in  a  sense 
that  will  be  defined.  In  Theorem  6.7  for  the  case  where  the  true  conditional  density  stays  close 
to  an  exponential  family  of  densities  we  show  that  the  error  of  the  estimate  given  by  projection 
particle  filtering  is  bounded. 

To  use  nonlinear  filtering  methods  for  carrier  phase  differential  GPS,  one  should  be  able  to 
include  the  integer  ambiguity  resolution  in  these  methods.  In  this  report  we  present  some  simulation 
results  which  show  that  particle  filtering,  with  minor  modifications,  is  capable  of  resolving  integer 
uncertainties  present  in  a  problem  similar  to  carrier  phase  differential  GPS.  One  problem  of  particle 
filtering  is  the  need  for  large  number  of  particles.  This  problem  is  even  more  important  for  the 
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cases  where  integer  uncertainty  is  present.  The  writers  are  optimistic  that  particle  filtering  for 
exponential  families  of  densities  is  more  suitable  for  nonlinear  filtering  with  integer  ambiguity. 

In  this  report,  Section  2  states  the  nonlinear  filtering  problem.  In  Section  3  we  review  the 
results  in  [11]  [12]  [13]  on  projection  filtering.  In  Section  4  we  explain  particle  filtering  and  we  state 
the  results  in  [8]  [9] .  In  Sections  5  and  6  we  introduce  a  new  particle  filtering  and  we  state  the  main 
results  of  this  report.  In  Section  7  we  apply  the  particle  filtering  to  a  nonlinear  system  with  integer 
uncertainty  and  we  present  the  simulation  results.  In  Section  8  we  discuss  the  future  research  on 
this  subject. 

2  Nonlinear  Filtering,  Problem  Setup 

Filtering  problems  consist  of  “estimating”  the  process  {x^}  (or  something  about  it)  given  the 
related  process,  {yt},  which  can  be  observed  [15].  The  observation  is  available  in  an  interval,  i.e., 
{ys,0  <  s  <  t}  and  the  function  of  the  state  is  estimated  at  time  t.  To  proceed,  we  need  to  give 
some  structure  to  the  concerned  processes. 

We  assume  that  all  stochastic  processes  are  defined  on  a  fixed  probability  space  (P,  F.  P),  and  a 
finite  time  interval,  [0,  T ],  on  which  there  is  defined  an  increasing  family  of  e-fields,  {Ft,Q  <  t  <  T}. 
It  is  assumed  that  each  process,  {xt},  is  adapted  to  Ft,  i.e.,  {x<}  is  ^-measurable  for  all  t.  We 
assume  that  {x^}  is  a  vector  diffusion  process  of  the  form 

x*  =  x0  +  /  f s(xs)ds+  f  Gs(xs)dws,  (1) 

Jo  Jo 

where  x*  £  7 Zn,  and  w t  £  lZq  is  a  vector  from  an  independent  Brownian  motion  process;  the 
second  integral  is  in  the  Ito  sense  [16],  and  the  function  f^(-)  and  the  matrix  Gt(-)  have  the  proper 
dimensions.  The  observation,  y^,  is  a  discrete  time  process  given  as  follows: 

ynT  =  hn(xnT)  +  vn,  (2) 

where  ynT  £  1Zd,  and  vn  £  1Zd  is  a  discrete  time  white  Gaussian  noise  process  with  zero  mean  and 
known  covariance  matrix.  The  state  dynamics  and  observation  equations  can  be  rewritten  formally 
as  follows: 

dxt  =  f)(xt)<if  +  Gi(xt)dwt,  given  the  distribution  of  xo 
y  nr  =  h„(xnT)+Vn 

The  noise  processes  {wt,  t  >  0},  and  {vn,  n  =  0,1,---}  ,  and  the  initial  condition  xo  are  as¬ 
sumed  to  be  independent.  We  use  Qt  and  Rn  for  the  covariance  matrices  of  the  processes  wt  and 
v„,,  respectively.  We  assume  that  Rn  is  invertible  for  all  n’s.  We  have  the  following  additional 
assumptions  [17]: 

A  2.1  [local  Lipschitz  continuity]  V  x,  x;  £  Br  and  t  £  [0,  T\,  where  Br  is  a  ball  of  radius  r,  we 
have 

llft(x)  -ft(x,)||  < 

|| Gt(x)QtGl(x.)  -  GtWQtGfWW  < 

A  2.2  [Non- Explosion]  There  exists  k  >  0  such  that 

xrft(x)  <  k(  1  + 
trace{Gt{yd)QtGj  (x))  <  k(  1 + 

M  t  £  [0,  T]  and  V  x  £  1Zn . 


&v||x  —  x7 1|,  and 

kr  1 1 X  —  X,||. 


(4) 


FI 

Ixll 


and 


(5) 
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Under  Assumptions  (A2.1)  and  (A2.2),  there  exists  a  unique  solution  {x<,  t  G  [0,  X]}  to  the 
state  equation,  and  x*  has  finite  moment  of  any  order  [17]. 

In  addition  to  these,  we  assume  that  the  probability  distribution  of  the  state  x^,  given  the 
observation  up  to  time  t.  irt(dx)  =  P(xt  G  rlxjy*),  where  y*  =  {yra,  i  =  1,  •  •  • ,  n,  nr  <  t},  has  a 
density  pt  with  respect  to  the  Lebesgue  measure  on  7 Zn.  Then  {pt,  t.  >  0}  satisfies  the  following 
partial  differential  equation  and  updating  equations  [12]: 

■§jpt  =  £'tPt  nr  <  t  <  (n+  l)r,  and 

Pnr  =  nPnT~ 

where 

4'(®)  =  -  Eti  AW®]  + 1  sfe-[“U]. 

[a^]  =  GtQtGj  > 

^n(x)  =  exp  (-|(ynT  -  hn(x))Ti?“1(ynr  -  h„(x)))  , 

and  cn  is  a  normalizing  factor. 

Except  for  the  linear  Gaussian  case,  and  some  very  special  nonlinear  cases,  solving  System  (6) 
constitutes  an  infinite  dimensional  filter  [10].  Therefore,  for  practical  problems  it  is  necessary  to  ap¬ 
proximate  the  conditional  density  in  (6).  In  the  next  section,  we  discuss  one  of  these  approximation 
methods. 


3  Projection  Filtering  on  Exponential  Families  of  Densities 

This  section  is  mainly  a  review  of  the  results  we  use  from  [12].  We  start  this  section  with  the 
definition  of  the  exponential  family  of  densities. 


Definition  3.1  Let  {ci,  •  •  •  ,cp}  be  affinely  independent  1  scalar  functions  defined  on  VJ' ,  and  as¬ 
sume  that  the  convex  set 

©o  =  I#  G  TZP  :  T (9)  =  log  J  exp  ^9Tc(x)'j  dx  <  oo|  , 
has  nonempty  interior.  Then, 

S  =  {pf,6),  0G0} 
p(x,9):=exp  6Tc(x)  —  T(6)  , 

where  0  C  @0  is  open,  is  called  an  exponential  family  of  probability  densities. 


We  denote  by  Sz  the  space  of  square  roots  of  the  densities  in  S  ,  i.e. ,  S?  =  {\/p(-,  0);  9  G  0}. 
If  p(-,  8)  G  S,  then  9)  G  L2.  The  functions  — p=  ^  dpgg^ ,  i  =  1,  •  •  •  ,p  form  a  basis  for  the 

tangent  vector  space  at  yjpf,  9)  to  the  space  S  2,  i.e.,  the  tangent  space  at  \/pf,  9)  is  given  by  [19]: 


L 


VpM) 


S 2  =  span 


1 

1 2  VpW) 


dpf,  8) 


dp{-,8) 


2y^M)  98 p 


(7) 


1{ci,  •  •  • ,  cp}  are  affinely  independent  if  for  distinct  points  xi,  X2,  •  •  •  ,  xp+i,  Ef-i  ^ic(xi)  =  0  and  =  0 

implies  Ai  =  A2  =  •  •  •  =  Ap+i  =  0  [18]. 
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The  inner  product  of  any  two  basis  elements  is  defined  as  follows 

1  dp(-,6)  l  dp(-,8)  \  _  l  C  1  dp(x,8)  dp(x,8) 

-  A  J  U’-X- 


2 y/p(-,8)  "i  ’  2v/p(-,0)  d63 


4  J  p(x,9)  ddi 

=  \9ij{0) 


dOi 


(8) 


It  can  be  easily  seen  that  g{6)  =  ( gijifi ))  =  (E[ciCj\  —  £'[cj]£,[cj])  is  the  Fisher  information  matrix 

of  p(-,0)- 

Any  member  of  L2  can  be  projected  to  the  tangent  space  L  jp^  ^<$2  according  to  the  following 
projection  formula 


n# :  L2  A  V 


dp(-fi) 


dOi 


dp(-,8) 


(9) 


'3  /  2v/^0)  90; 


Projection  filtering  seeks  a  solution  pt  for  (6)  that  lies  in  S.  Of  course,  this  solution  is  only  an 
exponential  density,  but  we  hope,  by  choosing  the  proper  family,  to  keep  the  approximation  error 
small  (in  the  L2  sense). 

If  we  consider  the  square  root  of  the  density  in  (6),  we  get 

dyfpt  _  1  dpt  _  1 


dt 


2y/Pt  9t  ‘Is/Pt 


(10) 


Define  at)e  =  We  assume  that  for  all  6  G  0  and  all  t  >  0,  Ep(.^{\atfl\2}  <  00,  which 

implies  that  is  a  vector  in  L‘>  for  all  0  G  O  and  all  f  >  0  [121. 

Now  assume  that  in  equation  (10),  for  {y/pt,  t  >  to},  starting  at  time  nr  from  the  initial 
condition,  yjpnr  =  \Jp(-,  0nr)  G  S2  for  some  6nT  G  0.  Under  these  assumptions,  the  right  hand 
side  of  (10)  is  in  L2,  which  can  be  projected  into  the  finite  dimensional  tangent  vector  space 
L  — -S? .  The  propagation  part  of  the  projection  filter  for  the  exponential  family,  S,  in  the 

interval  [nr,  (n  +  l)r),  is  defined  as  the  solution  to  the  following  differential  equation  in  the  same 
interval: 


d\/pt{-,  0t) 
dt 


=  n 


9t 


2  Vpt^A) 


(11) 


We  also  assume  that  h„(x)  in  equation  (2)  is  time  invariant,  i.e. ,  hn(x)  =  h(x),  and  the 


components  of  h(x),  hl (x) ,  and  ||h(x)  || are  linear  combinations  of  q(x),  i  =  1,  ■  ■  • ,  p: 
7jlih(x)ll|-i  =  Aic*(x)  and  /rfc(x)  =  ^  A^q(x),  k  =  1,  -  ■  ■ ,  rf 


(12) 


1=1 


i=  1 


where  ||x||^4  =  V xTdx.  Then,  if  vn  is  stationary  with  the  covariance  matrix  Rn  =  R,  the  likelihood 
function  \kn(n)  can  be  written  as  follows: 


^n(x)  = 


exp(-^(y nrR  1ynr))exp(-l(hT(x)ii  xh(x))  +  (ylTR  1h(x))) 

An  exp  (-  E  A?Ci(x)  +  E  (E  \^nr)Ci(X)  )  > 

V  i= 1  k=  1  i=l  J 


(13) 


where  znr  =  yJlTR~ 1 ,  and  An  is  a  constant  depending  on  ynT.  Therefore,  the  coefficient  \Hn(x)  is 
a  member  of  exponential  family  of  densities.  This  family  is  closed  under  multiplication.  Using  all 
of  these  facts,  we  can  present  the  following  theorem  [12]: 
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Theorem  3.2  [Brigo  1996]  For  system  (3),  where  w;  is  a  Brownian  motion  process  with  covari¬ 
ance  Qt  and  v*  is  a  white  Gaussian  noise  with  covariance  R,  we  assume  (A2A)  and  (A2.2)  to  be 

p  p 

true.  We  also  assume  that  ^||h(x)||^_1  =  A^Cj(x),  hk(x)  =  (x) ,  for  k  =  1  ,---,d,  and 

i= 1  i= 1 

Ep(.:e)  II  C,p,(l)f0})  l|2  <  oo,  ye  G  0,  yt  >  0.  Then  for  all  9  G  0,  and  all  t  >  0,  Tin  ~rf; =  is  a  vector 
on  the  exponential  manifold  S  2.  The  projection  filter  density,  p ^  =  Pt(-,9t)  is  described  by 


dy/pd-M 

dt 


Prar(*  ?  @nr)  — 


""•fiSS)’  nr<t<(n+  l)r 

cn  (y nr  )Pnr  ~  ('  ?  ^nr  “  )  ; 


and  t/ie  projection  filter  parameter  satisfies  the  following  combined  differential  and  stochastic  dif¬ 
ference  equations: 


g(9t)d9t  =  Egt{Ctc}dt,  nr<t<(n+l)r, 


9 nr  ~  9nr -  -  Ag  +  J2k= 1  ^0 zni 


where 


and  Ag  =  [A| ,  ■ 


1=1 


d 2 


2  ;  <)Xj()Xj  : 


,  Ap] 1 ,  i  =  0,  •  •  • ,  d,  and  zk  is  the  kth  component  of  z()T  =  R  1y„r. 


Henceforth,  we  shall  use  Eg  and  Ep^.  g^  9nr  and  9n,  and  pnT  and  pn,  interchangeably,  respectively. 

Remark:  The  differential  equation  for  9t  is  an  ordinary  differential  equation  with  the  vec¬ 
tor  field  g{9t)^1  Egt{Ctc}.  This  vector  field  should  be  computed  analytically.  If  the  analytical 
computation  of  this  vector  field  is  not  possible  an  off-line  numerical  computation  should  be  carried. 

As  can  be  seen  from  the  statement  of  the  theorem,  the  calculation  of  the  conditional  probability 
density  is  reduced  to  the  calculation  of  the  parameter  of  an  exponential  family.  But,  solving  the 
differential  equation  in  the  theorem  is  not  an  easy  task.  At  each  moment  g(9t)  and  Egt{Ctc} 
need  to  be  calculated.  This  imposes  a  heavy  computational  load.  In  this  report,  we  introduce  a 
Monte  Carlo  method  to  calculate  the  parameter  of  the  exponential  family  with  a  more  affordable 
computational  load. 

Although  projection  filtering  gives  a  better  solution  than  EKF,  there  is  no  known  error  bound 
with  which  we  can  compare  the  distance  between  the  real  density  and  the  density  given  by  the  pro¬ 
jection  filter.  In  the  next  section  we  review  particle  filtering  as  an  alternative  to  optimal  nonlinear 
filtering. 

Remark  :  The  assumption  on  h„. ( • )  and  Rn  in  this  are  made  only  to  ensure  that  is  in  the 

family  of  exponential  densities.  These  assumptions  can  be  relaxed  if  \kn(-)  is  guaranteed  to  stay  in 
the  family. 


4  Particle  Filtering 

Consider  either  the  continuous  dynamics  and  discrete  observation  in  (3)  or  the  discrete  case, 

xn+i  =  fn(xn)  +  Gn(xn)wn,  given  the  distribution  of  x0 
yn  =  hn(xn)  +  vn. 

We  assume  that  in  both  cases,  the  initial  distribution  for  xo  is  given.  The  propagation  of  the 
conditional  density,  at  least  conceptually,  can  be  calculated  as  follows  [10]: 
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•  Step  1  .  Initialization: 

p0(xo|yo)  =  p(x0). 

•  Step  2  .  Diffusion: 

p(n+i)_(xn+i|3^n)  =  J  p(xn+i  |xn)pn  (xre \yn)(htn, 

where  yn  =  {yi,y2,  •  •  •  ,yn}- 

•  Step  3  .  Bayes’  rule  update: 

p(yn4u|xn+i)p(n+i)_(xn+i|;yn) 

P(n+ 1)  Xn+1  ^n+1  J p(yn+1\xn+l)p^+i)_(xn+l\yn)dxn+l 


•  Step  4  .  n  <—  n  +  1;  go  to  Step  (2). 

The  conditional  density  given  by  the  above  steps  is  exact,  but  in  general  it  can  be  viewed  as  an 
infinite  dimensional  filter,  thus,  not  implement  able.  Particle  filtering,  in  brief,  is  an  approximation 
method  that  mimics  the  above  calculations  with  a  finite  number  of  operations  using  the  Monte 
Carlo  method.  The  procedure  for  particle  filtering  is  as  follows  [20,  8]: 

Algorithm  4.1  Particle  Filtering 

•  Step  1  .  Initialization 

o  Sample  xj,  •  ■  • ,  x^,  N  i.i.d.  random  vectors  with  the  initial  distribution  To(x)- 

•  Step  2  .  Diffusion 

o  Find  x*+1,  •  •  • ,  x^+1  from  the  given  x* ,  •  •  • ,  x^;  using  the  dynamic  rules: 

dxt  =  it(x.t)dt  +  Gtfx.t)dwt,  nr  <  t  <  (n  +  l)r 
or 


xri+l  —  fn(xn)  T  Gn(xn)vr 
Step  3  .  Find  the  empirical  distribution 


1 


N 


p<«+ i)-(x)  =  vS4i+1(x) 

3= 1 


•  Step  4  ■  Use  Bayes’  Rule 


•  Step  5  .  Resample 


P, 


N 


(™+i) 


(*)  = 


N 

jf  E  (X)  '  ^n+l(X) 

j  =  1  Xn+1 


1  N  , 

TV  E  ^ 

j=l  "+1 


(Xn+l)  •  ^n+l(x^+iJ 
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o  Sample  x*+1,  •  •  • ,  x^+1  according  to  -P^^+^x) 


•  Step  6  .  n  <—  n  +  1;  go  to  Step  (2). 

where  Sv(w)  =  1  if  w  =  v  and  0  otherwise,  and  \kn(x)  is  the  conditional  density  of  the  observation 
yn  given  the  state  x. 

It  is  customary  to  call  xT*( ,  •  •  ■ ,  x^  particles.  In  the  next  few  lines,  we  try  to  explain  in  words 
the  evolution  of  these  particles  using  the  above  algorithm. 

Let  x* ,  •  •  • ,  x^  be  the  distinct  particles  at  time  n  before  incorporating  the  observation  at 
time  n.  The  probability  of  each  particle  is  that  is,  is  uniformly  distributed.  After  using  the 
observations,  the  conditional  probability  of  each  particle  changes.  Some  will  have  small,  and  some 
large  probabilities.  Therefore,  in  the  process  of  resampling,  it  is  very  likely  that  some  particles  will 
never  be  used  and  instead  some  other  particles  (with  high  probabilities)  will  be  sampled  more  than 
once.  Therefore,  after  resampling,  some  particles  have  repeated  versions,  but  in  the  diffusion  phase 
they  go  through  different  paths  and  at  the  end  of  the  diffusion  phase,  it  is  very  likely,  we  would 
have  N  distinct  particles.  This  automatically  makes  the  approximation  one  of  better  resolution  in 
the  areas  where  the  probability  is  higher. 

In  [8]  it  is  proved  under  some  conditions  that 


lim  E 

A/"— XX) 


Epn(f(x)) 


=  0 


(15) 


for  every  bounded  Borel  test  function,  /(•). 

One  problem  in  using  the  particle  filtering  method  is  the  computational  cost.  In  particular,  for 
a  high  dimensional  system,  getting  reasonable  accuracy  means  using  a  large  N ,  which  results  in  a 
heavy  computational  cost.  In  the  next  section,  we  propose  a  method  that  can  reduce  the  number 
of  particles  for  a  certain  class  of  problems. 


5  Particle  Filtering  for  Exponential  Families  of  Densities 

In  the  previous  section,  we  saw  two  approximation  methods  for  nonlinear  filtering.  In  the  particle 
filtering  method,  we  saw  that  the  conditional  distribution  is  approximated  by  the  empirical  dis¬ 
tribution.  Unlike  the  empirical  distribution,  in  most  cases,  the  actual  conditional  distribution  is 
smooth.  Intuition  suggests  that  if  we  have  prior  knowledge  of  some  properties  of  the  distribution, 
we  can  improve  on  the  quality  of  the  estimates  over  just  using  the  empirical  distribution.  In  this 
section  first,  we  assume  that  the  conditional  density  lies  in  an  exponential  family  of  densities.  We 
will  see  that  with  this  assumption,  we  can  show  the  convergence  of  the  approximated  density  to 
the  actual  one.  Later,  we  relax  this  assumption  and  we  only  require  that  the  conditional  density 
stay  close  to  the  exponential  family  of  densities.  We  prove  that  the  error  of  the  estimate  for  the 
latter  case  is  bounded. 

For  System  (3),  we  assume  that  the  probability  density  of  x<,  given  the  observation,  is  in  a 
family  of  exponential  densities  S  2. 

With  this  assumption,  the  proposed  algorithm  is  as  follows: 

Algorithm  5.1  Particle  Filtering  for  an  Exponential  Family  of  Densities. 

2This  assumption  is  rather  strong.  We  will  drop  this  assumption  later,  and  we  will  only  assume  that  there  exists 
a  known  family  of  densities  that  approximates  the  real  density  well,  i.e.,  with  acceptable  accuracy. 


•  Step  1  .  Initialization 

o  Sample  Xg,  •  •  • ,  Xq ,  N  i.i.d.  random  vectors  with  the  density,  po(x). 

•  Step  2  .  Diffusion 

o  Find  x*+1,  •  •  • ,  x^+1  from  the  given  xf,  ■  ■  ■ ,  x^  ,  using  the  dynamic  rule: 


dxt  =  ff (xt)dt  +  Gt(xt)dwti  it  <t  <  (i  +  1  )r 
Step  3  .  Find  the  MLE  of  fn+1)-  given  x^+1,  •  •  • ,  x^+1  [21] 


N 

\n+ 1)-  =  argmax]Jexp(0rc(^+1)  -  T(0)) 

i= 1 


•  Step  4  ■  Use  Bayes’  Rule 


p{x,d{n+1)) 


eXP  +1)-C(x) 

-T(0(n+1)-) 

')  ®n+l(x) 

/exp(^+1)_c(x) 

-T(fn+1)-) 

)  ^n+l(x)c?X 

•  Step  5  .  Resample 

o  Sample  x^+1,  •  •  • ,  x^+1  according  to  p(x ,  6n+ 1). 

•  Step  6  .  n  <—  n  +  1;  go  to  Step  (2). 

To  generate  x^+1,  •  •  • ,  x^+1,  a  Gibbs  sampler  can  be  used  [22],  This  brings  an  extra  computational 
cost,  which  should  be  taken  into  account  when  choosing  Algorithm  5.1  over  Algorithm  4.1. 

It  is  instructive  to  discuss  the  structure  of  the  ML  estimator.  We  are  going  to  use  this  structure 
for  the  proof  of  convergence. 

Let  xln,  ■  ■  ■ ,  x][  be  the  value  of  the  particles  right  before  the  measurement  at  time  n.  The  MLE 
of  6n ,  6n,  satisfies  the  first  order  necessary  condition 

\'  Ar/xCi(X)eXP(^C(X))dx  n 

h'  "  /Xe*p(#c(x))d* 

Therefore,  we  get 


1 

for  j  =  !,•••, p  .  (16) 


Equation  (16)  says  that  the  sample  average  of  Cj(x)  and  its  probabilistic  average,  evaluated  at  6n , 
should  be  equal.  The  MLE  of  8  is  the  solution  to  the  system  of  equations  in  (16).  Let  Fj{9)  be  as 
follows: 


F^e) 


if  M  /  Cj(x)  exp((9rc(x))rfx 
Nfrf  n>  f  cxp(8Tc(x))dx 


j  =  !,•••,  P- 
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For  simplicity  we  drop  the  index  n  from  9n.  It  is  easy  to  see  that 

dF 

—q£  =  ^e(ci(x)cj(x))  -  Eo(ci(*))Ee(cj(x)). 

This  shows  that  (— =  g(9),  where  g(9)  is  the  Fisher  information  matrix  of  the  exponential 
density  at  9.  Since  c,;(x),  i  =  1  ,■■■  ,p  are  affinely  independent  g(9 )  >  0 ,\/9  £  0.  Therefore,  (16)  is 
the  necessary  and  sufficient  condition  for  optimality. 

In  the  next  few  pages,  we  prove  the  convergence  of  the  MLE  of  9n,  9n ,  to  9n  in  the  mean  square 
sense. 

In  each  iteration  the  proposed  algorithm  starts  from  the  density  Pq  (x*|y*),  t  =  rn,  where  9t 
is  the  best  estimate  9t  according  to  the  algorithm.  After  a  full  iteration,  the  algorithm  yields  9t+i 
which  is  the  best  estimate  of  6t+  \  ■  The  error  in  9t+ 1  is  a  combination  of  the  series  of  possible 
errors  for  which  we  want  to  find  an  upper  bound.  The  first  source  of  error  is  the  error  in  Of  .  which 
will  propagate  even  if  no  other  error  is  considered.  The  other  source  comes  from  the  fact  that  in 
each  iteration  new  particles  are  resampled  based  on  the  estimated  density  which  is  different  from 
the  actual  density.  Finally,  the  last  source  of  error  comes  from  the  discretization  of  the  stochastic 
dynamics  of  the  system.  We  want  to  emphasize  that  here  we  assume  \Iin(x)  =  exp(— ^(ynr  — 
hn(xr)r))ri?71(yrtT  —  hn(xnr)))  lies  in  the  chosen  family  of  densities.  Therefore,  no  other  error  is 
added  to  the  estimate  because  of  the  Bayes’  correction. 

We  recall  the  following  fact  [21]: 

Fact  5.2  For  the  family  of  densities  S  with  probability  density 

P(x,0)  =  exp(6»rc(x)  -  T(0)), 

the  Fisher  information  matrix  g (9)  =  (E(ci(x)cj(x))  —  E(cifx))E(cj(x)))itj  is  positive  definite.  Also 
the  log  likelihood  function 

1(9)  =  tfTC(x)  -T(0), 

is  strictly  concave.  Therefore,  for 

Cj(x)  =  Ee[cj{x)\,  j  =  l,---,p, 

if  a  solution  exists 3,  it  is  unique.  In  addition  if  x  i ,  •  •  ■ ,  x^r  are  N  i.i.d.  random  variables  distributed 
according  to  p(x,  9),  then  the  MLE  of  9,  9 n,  is  asymptotically  normal,  i.e. 

„  N 

9 N  =  arg  max  II  Kxo  0)  , 

_  e  i= l 

\/N(9n  —  9)  ~  M(t),g-1(9)). 

Using  this  fact,  it  is  easy  to  see  that 

E  9N  —  9  trace(g~1(9)), 

therefore,  when  N  — >  oo,  9n  — *  9  in  the  m.s.  sense.  On  the  other  hand,  9n  is  the  solution 
to  (16).  Using  the  strong  law  of  large  numbers  [23],  when  N  — >  oo  the  LHS  in  (16)  goes  to 

3In  [18]  it  is  shown  that  ii  N  >  p,  the  solution  exists  almost  surely. 
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Eg(cj(x)),  j  =  1  ,■■■  ,p,  with  probability  one.  In  other  words,  the  solution  to  (16)  when  the  LHS 
is  the  exact  Eg(cj(x)),  j  =  1  gives  the  exact  solution  for  6.  Using  this  argument,  one 

can  expect  that  by  finding  a  good  estimate  of  the  left  hand  side  of  (16),  a  good  estimate  of  0  is 
accessible.  In  each  iteration  of  the  algorithm  presented  in  this  section  the  estimate  of  the  LHS  of 
(16)  is  found  by  using  the  Monte  Carlo  method  and  the  approximate  solution  for  the  stochastic 
differential  equation  (3). 

To  approximate  the  solution  to  the  stochastic  differential  equation  (3),  we  employ  the  method 
used  in  [24],  In  the  following,  we  review  this  method  briefly.  The  stochastic  differential  equation 
in  (3)  can  be  rewritten  as  follows: 


dxf  =  it  (xt)  dt  +  ^2  grt  {xt)  dWf ,  (17) 

r— 1 

where  g[(-)  is  the  rth  column  of  the  matrix  Gt(- ),  and  is  the  rth  component  of  W(.  We  introduce 
the  operators 


A  ru 


u, 


Lu  = 


nEEE*: 


o'- 


r= 1  i= 1  j=  1 


J  <9xj<9xj 


U, 


where  (a  ,  ^ )  =  E  aiJET-  Then,  the  approximate  solution  for  the  stochastic  differential  equation 
'  '  i= 1  1 

can  be  written  as  follows: 


xfc+1  =  xfc  +  E  gf&h*  +  f tkh  +  E  E  (Argr)tfe  ikh+ 

r= 1  r=l  i=l 

hi(Lsr  +  Kf)tJlhi+(L() 

r=  1 


(18) 


where  h  is  the  step  size  and  the  coefficients  g[fc  ,  f tk,  (Aig *•)*,,  etc.,  are  computed  at  the  point 
(tfc,Xfc).  Also,  the  sets  of  random  variables  ££,  ££  are  independent  for  distinct  k  and  can,  for  each 
k,  be  modeled  as  follows: 


«« =  \ee  -  \njCV, 


lij 


-1  ,  'i  <  J 

1  ,i>j 


and  £*  and  ^  are  independent  random  variables  satisfying 

Eti  =  E(f  =  E$  =  0,  E%  =  1,  Eif  =  3, 

E(j  =  ECf  =  0,  E(]  =  Cf  =  1. 


In  particular,  Ci  can  be  modeled  by  the  law  P  (C  =  0)  =  |,  P  ^  =  \/Sj  =  g,  and  P  \  £  =  —  V3 
and  (j  can  be  modeled  by  P  (C  =  —1)  =  P  (C  =  1)  =  |- 
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Definition  5.3  We  say  that  a  function  u(-)  belongs  to  the  class  F ,  written  as  u  E  F ,  if  we  can 
find,  constants,  k  >  0,  and  k  >  0,  such  that  for  all  x  E  lZn ,  the  following  inequality  holds: 

Hx)ii  <  jfe(iH-nxir). 


Before  we  present  our  results  we  need  to  specify  the  probability  space  in  which  the  random 
variables  are  defined.  As  we  mentioned  before,  the  stochastic  process  associated  to  the  dynamics 
and  the  observation  equation  are  defined  on  a  fixed  probability  space  (Q,F,P),  the  expectation 
associated  to  this  probability  space  is  denoted  by  E.  In  Algorithm  5.1  the  generated  particles 
form  a  Markov  process.  Similar  to  section  2.2  of  [8]  we  denote  the  probability  space  associated  to 
this  process  by  The  subindex  y  is  used  to  emphasize  that  the  probability  measure 

is  conditioned  on  the  observation  y.  Another  set  of  random  variables,  £*,£*,  are  defined  for  the 
approximation  of  the  stochastic  differential  equation  (17).  We  denote  the  probability  space  asso¬ 
ciated  to  these  random  variables  by  (Q," ,  F" ,  P") .  The  expectation  associated  to  this  process  is 
denoted  by  E" .  Finally  we  define  (Cl,F,P),  where  Cl  =  Q  x  Q1  x  Q,11  and  F  =  F  x  F'  x  F" .  For 
every  Cj  E  fl  we  define  u)  =  (cu,  u/,  w"),  then  for  every  A  E  F,  B  E  F' ,  and  C  E  F"  we  define  the 
probability  measure  P(A  x  BxC)  =  fA  (  fc  P^(B)dP" dP(uj).  The  expectation  with  respect 

the  probability  measure  P  is  denoted  by  E. 

The  following  theorem  summarizes  the  weak  approximation  results  for  (18). 

Theorem  5.4  [  Milstein  [24]]  Suppose  (A2.1)  from  Section  (2),  and  suppose  that  the  functions 
f(0,  gr('),  r  =  1,  -  ■  ■ ,  q  together  with  the  partial  derivatives  of  sufficiently  high  order,  belong  to  class 
F .  Also,  suppose  that  the  functions  A*gr,  Lgr,  Arf,  and  Lf  grow  at  most  as  a  linear  function 
in  |jx||.  Then,  if  the  function  u(-)  and  all  its  derivatives  up  to  order  6  belong  to  class  F ,  the 
approximation  (18)  has  the  order  of  accuracy  2,  in  the  sense  of  weak  approximation,  i.e., 

||-S«(xo,x0  (4))  -  Eu  (x0ix0  (4))  II  <  A7?2,  tk  E  [0,T], 


where  K  is  a  constant  and  xq)Xo(-)  and  xoiXo(")  are  the  exact  and  approximate  solutions  for  the 
stochastic  differential  equation,  respectively. 

The  Monte  Carlo  approximation  of  Eu(x o.X(,  (4))  brings  another  error  term.  The  combination  of 
these  errors  can  be  expressed  as  follows: 


Eu  (x0,xo  (4)) 


Eu(x0,x o  (4)) 


Eu  (x0,xo  (4)) 


+ 


Eu  (x0,x0  (4)) 


If  the  variance  of  u  (xo,Xo  (4))  is  bounded,  we  have 


E 


Eu  (xq.xo  (4)) 


1 

N 


N 

5> 

i= 1 


(xo,xj  (4)) 


<  Kh2  + 


k' 


(19) 


where  K  and  k'  are  constants,  and  h  is  the  step  size  for  the  approximation  of  the  solution  of  the 
stochastic  differential  equation. 

The  next  lemma  relates  the  approximate  solution  to  the  stochastic  differential  equation  and  the 
estimate  of  the  parameter  6.  This  lemma  is  the  main  building  block  for  our  result  in  this  section. 
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Lemma  5.5  For  the  stochastic  differential  equation 

dxt  =  it(xt)dt  + Gt{xt)dwt,  x0,  te[0,tf], 

assume  that  f* (-) ,  Gt(-)  are  such  that  for  the  Brownian  motion,  wj,  the  probability  density  of  the 
state  x*  lies  in  the  family  S  for  0  bounded,  with  g{9)  >  01  for  some  d  >  0.  We  also  assume  the 
conditions  in  Fact  5.2  and  in  Theorem  5.4  with  c(x)  replacing  it(x).  Then,  there  exist  k\  and  k,2 
such  that 

E[\\9t-9t\\\<k1h2  +  t€[0,tf]  (20) 

where  9t  is  the  estimate  of  9t,  and  N  and  h  are  the  number  of  particles  and  the  time  step,  respec¬ 
tively. 


Proof:  Let  9 o  be  the  initial  condition  for  9.  At  t  =  0,  N  independent  initial  conditions  are 
generated  based  on  the  density  p(x,  9q),  and  the  approximation  method  (18)  is  applied.  From  (19) 
we  know  that: 

1  N  k' 

E\\Eetc  (xt)  -  £  c  (* 0  II  ^  Kh2  +  ^172  • 

i=  1 

On  the  other  hand,  from  (16),  we  know  that  9  is  a  solution  to  the  system  of  equations 

1  N 

^  £cj(xt)  =  E^Cjfxt)),  for  j  =  1,  •  •  •  ,p. 


N 

From  Fact  5.2,  the  solution  is  exact  if  we  replace  jr  Cj(xJ)  by  Egt(cjfx.t)).  Subtracting  the 

i= 1 

term  Egt(cjfx))  from  both  sides  of  the  above  equations  and  using  the  vector  form  for  it,  we  get 

1  N 

jy  £  C(^i)  “  Eet{c(xt))  =  Egt( c(xt))  -  E0t(c(xt)). 

i= 1 

On  the  other  hand,  we  know  that  Eg( c(x))  is  a  differentiable  and  one  to  one  function  of  9  (  see 
Fact  5.2).  The  derivative  of  this  function,  g(9),  is  positive  definite  and  by  assumption  g{6)  >  dF 
Therefore,  3a  >  0  such  that 

\\0t-dt\\  <  a\\Egt(c(xt))  -  E§t(c(xt))\\ 

N 

=  a\\Egt(c(xt))  -  E  c(xj)||. 

i= 1 

Taking  the  expectation  on  both  sides  of  the  inequality  we  have 

E\\et-et\\  <  a  E\\E  £  c(x|)  -  Egt(c(xt))\\ 

i=  1 

<  a  ^ KK 2  + 

=  fa  h2  +  ^  o 

Now  we  are  ready  to  present  the  main  result  of  this  section. 
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Theorem  5.6  For  System  (3)  assume  that  f [(•),  Gt(-),  and  h(-)  are  such  that  for  the  Brownian  mo¬ 
tion  wt,  and  the  Gaussian  noise  vn,  the  conditional  probability  density  of  the  state  xt;  conditioned 
on  the  observations,  lies  in  the  family  S  for  0  bounded  and  for  t  £  [0,  T\.  Also  assume  the  condi¬ 
tions  in  Fact  5.2  and  in  Theorem  5.4  with  c(x)  replacing  ufx).  Then,  if  g^1  {Of)  Egt  (£tc(x))  is 
Lipschitz  with  Lipschitz  constant  L  and  g(0)  >  •&!,  there  exist  l\  and  I2  such  that 

n-l  /  7  \ 

E\\0n  -  0n ||  <  Y  exp (Lir)  y>ih2  +  ,  nr  £  [0,  T], 

where  0n  is  the  estimate  of  0n,  and  N  and  h  are  the  number  of  particles  and  the  time  step,  respec¬ 
tively.  This  inequality  implies  convergence  of  the  estimated  parameter,  0n,  to  the  true  parameter, 
0n,  as  h  — >  0  and  N  — >  00. 

Proof:  Let  0t  and  Of  be  the  actual  and  the  estimated  values  of  the  parameter  of  the  density  at 
time  t  =  nr,  respectively.  At  time  t  =  (n+  l)r  the  error  in  the  estimate  of  Op  is  a  combination  of 
the  error  in  the  estimate  in  Of  and  the  error  added  in  the  time  interval  [t,t]. 

If  the  conditional  density  stays  in  the  exponential  family  of  densities,  0t  has  to  satisfy  the 
following  differential  equation: 

0  =  g-1  {0)  Egt  {Ct.c  (x))  dt ,  nr  <  t  <  (n  +  1)  t. 

Let  O.1  be  the  estimate  of  Op ,  if  the  error  due  to  resampling  and  the  approximation  of  the 
stochastic  differential  equation  solution  is  not  taken  into  account  in  the  interval  [t,  t  ]  (i.e.  Of  is 
computed  from  the  above  ordinary  differential  equation  starting  at  Ot),  then 

II °t'  ~  Ml  ^  II 6t’  -  Ml  +  II -  ^'11- 

By  the  assumption  of  the  theorem,  g~l  ( 0 )  Egt  (Ct c(x))  is  Lipschitz  with  Lipschitz  constant  L , 
then  by  continuity  of  the  solution  of  the  differential  equation  with  respect  to  the  initial  condition 
[25],  we  know  that 

Of  -  Of  <  0t  -  0t  e , 

therefore, 

E  Of -Of  <E  0t-0t  eL(t'_t). 

Also  from  the  Lemma  5.5,  3ki{t)  and  k2(t  )  such  that 

i#,' -  Ml]  <  M*>2  + 

therefore, 

E\\0p  -0t'\\  <E\0t-0t\eL^-^ +  k1{t)h2+1^Jl. 

The  observation  noise  v„  and  the  function  h(-)  are  such  that  Bayes’  Rule  does  not  introduce  any 
further  error  in  the  estimate  of  Op.  More  precisely,  'Lre(x)  is  assumed  to  be  a  member  of  S.  This 
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implies  that  after  applying  Bayes’  Rule  to  p(x,  6t>)  and  p(x,  9tr)  parameters  9ti  and  9t/  are  shifted 
with  the  same  vector  and  therefore,  || 9t+>  —  9t+/ 1|  =  || 9j  —  9 j  |.  Here  t,+  represents  the  time  right 
after  Bayes’  correction.  Therefore,  starting  from  the  initial  condition  9q  we  get 


where 


E\\9n 


n—  1 

9 n  1 1  <  ^  exp(Lir) 
i= o 


nr  G  [0,  T] 


li  =  max ki (nr),  nr  G  [0,71,  i  =  1,2. 

n 


O 


Here,  we  would  like  to  make  a  few  remarks: 

•  The  result  of  Theorem  5.6  can  be  easily  extended  to  convergence  in  the  mean  square  sense. 

•  The  assumption  that  the  probability  density  stays  in  the  family  of  densities,  S,  does  not  seem 
very  realistic.  But  with  our  approach,  we  should  be  able  to  get  the  result  in  [12].  In  fact,  in 
[12]  the  evolution  of  the  density  is  forced  to  stay  in  the  family  at  every  single  moment.  In 
our  method,  we  only  force  the  density  to  be  in  the  family  at  the  end  of  each  full  iteration,  i.e. 
observation  epoch.  This  allows  the  estimated  density  to  be  closer  to  the  actual  density. 

•  In  [12]  the  observation  equation  is  considered  to  be  time  invariant.  Here,  the  time-varying 
nature  of  h„  (x)  does  not  complicate  the  algorithm.  It  surely  affects  the  assumption  that  the 
density  stays  in  the  family,  but  as  we  explained  earlier,  this  assumption  is  not  realistic  to 
begin  with,  and  it  will  be  dropped. 

•  If  u(-)  is  in  J7,  then 


lim  E\\E$uCx)  —  Eq*uCx)\\  =  0. 

N - >oo 

h - >0 

This  is  a  criterion  similar  to  the  one  used  in  [8] . 

6  Projection  Particle  Filtering  for  Exponential  Families  of  Densi¬ 
ties 

In  this  section,  we  drop  the  assumption  that  the  conditional  density  of  the  state  given  the  obser¬ 
vation  (6)  lies  in  the  exponential  family  of  densities,  S.  Also,  we  do  not  require  that  \kn(x)  is  a 
member  of  S.  Instead  we  make  other  assumptions.  First  we  need  the  following  definition: 

Definition  6.1  We  say  that  a  function  u(-)  belongs  to  the  class  written  as  u  G  J~kK,  for  fixed 
k  >  0  and  k  >  0,  such  that  for  all  x  G  lZn ,  the  following  inequality  holds: 

Hx)||  <&(l  +  ||x||K). 

The  next  two  assumptions  are  to  guarantee  the  existence  of  an  exponential  density  close  to  the 
true  conditional  density. 
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A  6.2  For  the  density  in  (6)  there  exists  an  exponential  family  of  densities  S  such  that  Vi  G 
[0,  T\,  Vu  G  tFkn  30)*  €  0*  and  e  >  0  such  that 

E\\EPt(u(x))  -  Fg*(n(x))||  <  e  ,  (21) 

where  0*  is  convex  4  and  compact. 

A  6.3  For  0*_  in  (A6.2)  and  Vl/n(x),  3T*(x)  such  that 

=  P(x,g;_)^(x) 

f  p(x,0*_)V*(x)dx 

is  in  the  family  S  for  some  0  G  0*  and  we  have: 

•  V0  G  0*  and  Vtt(-)  €  3e  >  0  such  that 

£  Ee^n(x)u(x)  _  Fg'h*  (x)n(x)  < 

Fg\Ifn(x)  Fg^*(x) 

•  Vn(-)  G  Tkn,  3e  >  0  such  that 

~  Ee*  ^*(x)u(x)  Fp  ^n(x)u(x) 

Fll  —  n _  _  n _ 1 1  ^  ^ 

Ee*  ^n(x)  Fp  ^n(x) 

n—  n 

From  Assumption  (A6.3)  it  is  clear  that  if  \H*(-)  satisfies  the  requirements  of  the  assumption  then 
c\f,*(-)  satisfies  the  same  requirements,  where  c  is  a  positive  constant.  Therefore,  without  loss  of 
generality  we  assume  that  \If*(-)  =  exp(aTc(-))  for  some  a  G  1ZP .  Using  Assumption  (A6.2),  we 
can  state  the  following  fact. 

Fact  6.4  V$i,$2  G  0*  andMu  G  T^k,  3F'i,F'2  positive  such  that 

||Fg1n(x)-Fg2n(x)||  <  K^d i-02||  (22) 


||di-d2 1|  <  F'2||Fg1c(x)  —  Fg2c(x)[ 

Proof:  To  prove  (22),  define  /u(0)  =  Egu{x )  for  u(-)  G  F/CK ■  Then 

d 


(23) 


dO 


■/«(0)  =  Egd(x)u(x)  -  EgCi{x)Egu(x). 


4It  is  easy  to  see  that  the  assumption  of  convexity  is  very  natural.  Assume  61,62  £  0*  then  f  exp(df  c(x))dx  <  00 
for  i  =  1,2.  Therefore,  using  the  Holder  inequality  we  have 


J  exp((ad^'  +  (1  —  a)0^)c(x))dx  —  f  (exp(0^c(x)))“  (exp((9^c(x))) U  dx 

<  ^/((ex P(efc(x)))  )  '  dxj  ^(exp(S|’c(x))) 

=  (J*  exp(0^c (x))dxj  exp(0^c (x))dxj 


l/l-o 


where  0  <  a  <  1. 
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Since  ||tt(x)||  <  k(  1  +  ||x||K)  and  9  G  0*,  where  0*  is  compact,  then  there  exists  a  constant  A  such 
that 

\\^4P-\\<A  \/u{-)eEkK  and  V0G0*. 
du 

Since  0*  is  convex  and  compact,  it  is  clear  that  3K\  independent  of  u(-)  such  that  /u(x)  is  Lipschitz 
over  0*  with  the  Lipschitz  constant  K\  [25]. 

Inequality  (23)  follows  from  the  fact  that  0*  is  compact  and  the  Fisher  information  matrix 
g{9)  >  dl  over  0*. 


o 


Denote  the  interior  of  the  set  0*  by  0(n^.  For  0^  we  can  state  the  following  fact. 

Fact  6.5  The  set 

A  =  ja  :  J  exp(aTc(x))  exp(0Tc(x))dx  <  oo ,V0  G  0*f^  and  a  G  FTP 

is  closed. 

Proof:  Assume  A  is  not  closed.  Therefore,  there  exists  a  converging  sequence  {a*}  C  A  with 
converging  point  a  ^  A,  then  39  G  0?^  such  that 

exp(aTc(x))  exp(6,Tc(x))dx  >  M,  VM  G  7 Z. 


Since  0-^  is  an  open  set,  3e  >  0  such  that  A fe(9)  G  0-^.  Also,  since  {«*}  is  a  converging  sequence, 
3k  >  0  such  that  ak  G  Afja).  This  implies  that  9\  G  0-^  where  9\  =  9  +  a  —  ak.  Therefore, 

J  exp(a^c(x))  exp(0fc(x))dx  <  oo. 

On  the  other  hand,  we  know  that 

exp(a^c(x))  exp(#fc(x))  =  exp(aTc(x))  exp(0Tc(x))  . 

This  is  a  contradiction,  therefore,  A  is  closed. 


o 


The  following  lemma  is  one  of  the  building  blocks  of  the  results  of  this  section. 

Lemma  6.6  For  9*_  and  VH*  (x)  defined  in  (A6.3),  andVu(-)  G  Tkn,  3  positive  numbers  £q,  k2,  Aq,  Aq 
independent  of  0*_  and  Vl/*(x),  such  that  39\.  92  G  0*  the  following  are  true. 

(a)  k1<\\EeK(x)\\<k2  V0  G  0*. 

(b)  )|F^(x)u(x)||<A;3  V0  G  0*. 

(c)  \\E9l ®*(x)«(x)  -  F02^(x)«(x)||  <  k4\\9i  -  92\\. 
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Proof:  Let  At  be  a  set  defined  as  follows 


M  =  {m  :  m  =  9\  —  62,  V6b,#2  £  ©*}• 

We  claim  that  M  is  compact.  To  prove  this  claim,  assume  {m*}  to  be  a  sequence  in  At,  i.e  m*  6  At. 
Also  we  assume  that  lim  m;  =  m.  We  know  that  there  exist  sequences  and  {$2,1}  such  that 

i — >00 

m,  =  Oil  —  62,1  and  0i,i,  02,i  £  0*.  Since  0*  is  compact  there  exist  converging  subsequences 
and  {62 ,i}  in  0*.  This  implies  that  m  =  9\  —  9 2,  where  6\  and  9 2  are  the  limits  of  the  subsequences 
{9 \ti}  and  {#2,*}-  Since  9\  and  6*2  £  0*,  then  m  £  At,  therefore  M.  is  closed.  Since  0*  is  bounded, 
At  is  bounded  and  therefore,  it  is  compact. 

We  define  set  A\  as  follows: 

A\  =  ja  :  J  exp(aTc(x))  exp(0Tc(x))dx  <  00  ,V0  G  0*  and  a£lZp 

It  is  clear  that  A\  C  A.  As  we  mentioned  before,  without  loss  of  generality  we  can  assume 
dp^x)  =  exp(aTc(x))  and  from  Assumption  (A6.3)  it  is  clear  that  a  £  Af ]  AL  Since  AC |  AT  and 
0*  are  compact  we  have 

min  min  llE'e'I'*  (x)||  <  HFL’I'*  (x)  ||  <  max  max  \\Eg^>*  (x)||. 

dee*a£Af)M 

In  other  words  (a)  is  true  with  k\  =  min  min  \\Eg^>*  (x)  II  and 

0e©*  Q£.4  p|  m 

ko  =  max  max  (x)||.  Similarly,  since  u(-)  £  EkK,  (b)  is  true. 

0e©*  aeAf]M 

Using  the  above  argument  and  the  argument  in  Fact  6.4,  we  can  show  that  ||^U04'*(x)u(x)|| 
is  uniformly  bounded  and  since  0*  is  convex  and  compact,  then  (c)  is  true  [25]. 

o 

In  the  following  we  go  through  the  proof  of  the  theorem  that  we  state  later  precisely.  Assume 
9n  is  calculated  according  to  Algorithm  5.1  and  assume  p(x,  9n)  is  such  that  \/u  G  E^k 

E\\E$nu(x)  ~  e9*u(x)\\  <  5  ,  (24) 

where  (see  Assumption  (A6.2))  satisfies 

E\\EPnu(x.)  -  Ee*u(-x)\\  <  e.  (25) 

Using  the  density  p(x,  9n),  new  particles  x*,  •  •  • ,  x(^  are  generated.  The  approximate  solution  for 
the  stochastic  differential  equation  at  time  (n  +  1  )r  maps  these  particles  to  x*+1,  •  •  • ,  x.()r+  j .  From 
these  new  particles  9n+\  is  calculated.  From  (24)  and  (25)  we  have 

E\\EPnu{-x)  -  E§nu(-K)\\  <  d  +  e.  (26) 

We  define  the  function  r(x)  as  follows: 

r(x)  =  E"c(xn>x((n+  1))) 

where  x„,x((n  +  l)r)  is  the  approximate  solution  of  stochastic  differential  equation  (17)  at  time 
(n+  1  )r  with  the  given  initial  condition  x  at  time  nr  using  the  method  in  (18).  Since  according  to 
our  assumption  c  G  then  by  using  lemma  9.1  in  [24],  we  have 

||r(x)||  <  tf3(l  +  Hxf), 
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where  /I3  and  /+  only  depend  on  the  function  c(-)  and  the  dimension  of  x.  We  assume  that  r  £  Ekn. 
If  the  argument  of  r(-)  is  a  random  variable,  then  using  (26)  we  have 


E\\EPnv(-K)  -  %r(x)||  <  5  +  e. 


More  explicitly, 

E\\EPnE"[c(xn,x((n  +  l)T))\  -  EgnE”[c(-kn^((n+  l)r)]||  <  8  +  e. 

From  Theorem  5.4  we  have 


E\\EPnc(xn^((n  +  l)r))  -  EPnE"  c{±n^{(n  +  1)t))\\  <  I<4K2, 


for  some  K 4  >  0. 

Using  the  Monte  Carlo  method  to  calculate  the  i?Pnc(xn)X((n  +  l)r))  brings  another 
term  that  is  due  to  the  finite  number  of  particles  as  the  initial  conditions  for  method  (18). 
expectation  of  this  error  is  bounded,  i.e.  37+  >  0  s.t. 


N 


E\\ EenE”  c^nA{n  +  1)t))  ~  c(*n+((re  +  !)r 


< 


i=  1 


ih 


where  x'  are  distributed  according  to  p(x,9n).  Combining  (28),  (29),  and  (30)  we  get 
^l|£pnc(xn,x((n+  l)r))  -  jf  J2iLl  c(xn,x*((n  +  l)r))ll  < 


<5  +  e  +  IUh2  + 


Based  on  (A6.2),  we  know  that  3 0(n+1)-  such  that 


E\\E., 


P(n+ 1) 


_c (x)-E0»  c(x)  ||  <  e. 

(n+1) 


We  know  that,  if  x  (initial  condition  at  time  nr)  is  distributed  according  to  pn(x),  then  Ep 
£'Pnc(xrtiX((n  +  l)r)),  therefore,  from  (31)  and  (32)  we  get 


E\\E0*  ,  c(x) 

(n+1) 


1  K 

T7  c(xn,xi((n  +  l)r))ll  <  8  +  2e  +  K4h2  H  -p 
i=i  iV2 


Then  +n+i)-  given  by  Algorithm  5.1  satisfies  the  following  inequality 


E\\Ee*  c(x)  -  En  c(x)||  <  5  +  2e  +  K4h2  + 

(«+i)-  >+1)- 


#5 

IVs’ 


From  (A6.3)  we  know  that  3 6  £  0*  such  that 


Ee*  ^;+i(x)u(x)  Ep  tf„+i(x)u(x) 

(n+1)  (n+1) 

=  E 

Eeu(x)  -  Ep^+i)u(x) 

Ee*  ^+iW  Ep  'Um-i(x) 

(n+1)'  (n+1) 

;!  e. 


Since  6  satisfies  the  inequality  in  (A6.2),  we  can  choose  $*n+:L)  to  be  6 , 


i.e. 


(27) 

(28) 

(29) 

error 

The 

(30) 

(31) 

(32) 
c(x)  = 

(33) 

(34) 
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On  the  other  hand  we  have 


therefore, 


Eg*  ii(x)  —  Ea  u(x) 

%+l)  V  )  6  (n  +  l)  V  ' 


< 


Eg*  ^*+i(x)m(x) 

(n+l)- 

Efj  4>n+l(x)u(x) 

(n+l)” 

Eg*  *;+1(x) 

(n+l) 

Eg  4/n  +  l  (x) 

(n+1)- 

Eg*  «r;+1(x)w(x) 

(n+l)- 

Eg*  ^*+1(x)u(x) 

(n+l)' 

Eg*  K+lW 

(n+l) 

+n-^+l(x) 

(n+l) 

Eg*  ^+i(x)«(x) 

(n  +  l)' 

Eg  4>;+1(x)U(x) 

(n+l) 

Ee(  ,,_^+iW 

(n+l) 

+n_^+iW 

(n+l) 

^  *;+1(x)«(x) 

(n+l) 

4>n+i(x)u(x) 

(n+l)” 

Ee.  +1.+;+iW 

(n+l) 

Eg  4/n  +  l  (x) 

(n+l)' 

+ 


+ 


Eg*  ll(x)  —  Ex  u(x) 

%  +  l)  K  ’  °(n+l)  y  ' 


< 


II  Eg*  *;+1(x)u(x)|| 

W7*  CwR  *++d 

(n+l)-  ("+!)- 


E0*  ,,/n+lW-^  ^n+l(x) 

(n+l)  (n+l)- 


+ 


WEe, 

(n+l) 


Eg*  ®n+1(x)«(x)  -  E§ 

(n+l)  T  D(n+1)“ 


®n+l(XMX) 


+ 


Eg  *n+1(XMX) 

(n+l) 

Ex  'I'n+i(x)u(x) 

(n+l)" 

Ex  41*  ,  .  (x) 

(n+l) 

Eg  ^n+l(x) 

(n+l)" 

Using  Lemma  6.6  and  (A6.3)  we  get 

E\\E0L^Ax)  ~Eelr,+uUWW  -  k3k4^klk4E\\ flfn+1)_ 


(n+l)  v  '  y(n+ 1) 

Therefore,  from  (34)  and  Fact  6.4  we  get 


\n+ 1)-  II  +  e- 


£||0(„+i)-  -  0(n+ 1)-  ||  <  ^2  U  +  2e  +  K4hz  +  -f  • 

iV  2 


This  implies  that,  3ti,  (-2U3U4  >  0  such  that 


U||£+  m(x)  -  Ex  u(x)  II  <  ii<5  +  t2C  +  +  (-4^  2  . 


(n+l) 


(n+l) 


The  next  theorem  summarizes  our  result  in  this  section. 


Theorem  6.7  For  the  system  (3)  assume  (A2.1),  (A2.2),  (A6.2),  and  (A6.3).  We  also  assume 
the  conditions  in  Fact  5.2  and  in  Theorem  5.4  with  c(x)  replacing  u(x),  and  we  assume  r  e  Ekn- 
Then  in  Algorithm  5.1  with  approximation  (18),  if\/u(-)  €  EkK 

E\\E§nu(x)  ~  Eo*nu(*)\\  <  5 

then 
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E\\Ee*  u(x)  -  En  u(x.)\\  <  LiS  +  i2e  +  i3h2  +  i±N  2, 

(n+1;  w(n+ 1) 

for  some  l\,  l2,  l3,  la  >  0. 

In  Theorem  6.7  only  one  step  of  Algorithm  5.1  is  considered;  it  is  straightforward  to  then  use 
Theorem  6.7  repeatedly  for  the  time  interval  [0,  T],  where  T  =  Mr.  In  that  case,  || Eqqu(x)  — 
Eg*u(x)  ||  <  e>o,  then  204,0:2,0:3,  and  0:4  positive  such  that 

E\\E0*u(x.((ri)T))  -  £'^u(x((n)r))||  <  a\5o  +  a2e  +  a3h2  +  a4lV_1/2, 


for  0  <  n  <  M. 

7  Particle  Filtering  for  Nonlinear  Systems 
with  Constant  Integer  Uncertainty 

Consider  the  following  nonlinear  dynamics  and  observation 

dxt  =  it  (xt )  dt  +  Gt  (xt )  dwt 
y  nT  =  h  „(x(nr))  +  Jn  z  +  vn 

where  the  assumptions  and  the  dimensions  for  x^,  ynr  ,  W(.  and  vn  are  the  same  as  in  the  previous 
sections.  We  assume  that  z  is  a  random  integer  vector,  i.e.  z  €  Zm  and  Jn  has  the  proper 
dimension.  Vector  z  is  assumed  to  be  constant  in  time.  This  problem  can  be  set  up  in  discrete 
time  as  well.  In  this  case, the  system  dynamics  and  the  observation  can  be  written  as  follows: 

Xn+1  =  fn(Xn)  +  Gn(xn)wn 
yn  =  hn(xn)  +  Jnz  +  vn 

In  both  setups  we  assume  that  the  integer  uncertainty  affects  only  some  components  of  the 
observation,  and  other  components  are  unaffected  by  z.  The  affected  components  have  associated 
noise  components  in  vn  that  have  considerably  lower  energy.  In  other  words,  the  uncertain  compo¬ 
nents  of  ynr  (or  equivalently  yn)  would  be  considerably  more  accurate  than  the  other  components, 
if  the  integer  ambiguities  were  known.  This  suggests  that  an  accurate  estimation  of  z  can  increase 
the  accuracy  of  the  estimate  of  the  state  of  the  system  significantly.  With  this  explanation,  our 
treatment  of  z  is  clear.  From  the  state  dynamics  and  the  observation  equation  we  first  estimate 
z  and  then,  with  fixed  z,  we  use  regular  nonlinear  filtering  methods  to  estimate  the  state  of  the 
system  xt. 

We  augment  the  state  x^  with  the  integer  ambiguity  z.  Having  done  that,  the  state  dynamics 
and  the  observation  have  the  following  form: 


z  t 

Ynr 


Uxi) 

0 


dt  + 


Gt{xt) 

0 


dwt 


(35) 


=  hn(x(rar)  +  Jnz(rar))  +  vr 


We  assume  that  the  initial  distribution  of  (xq  ,Zq)t  is  known.  Now  the  state  dynamics  and  the 
observation  have  the  same  form  as  was  studied  in  Section  (4).  Therefore,  we  can  apply  particle 
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filtering  to  find  the  conditional  probability  distribution  of  the  augmented  state.  This  setup  is  a 
special  case  of  the  setup  in  Section  (4).  In  (35)  there  is  no  state  transition  for  zt,  therefore,  using 
particle  filtering  in  its  original  form  may  not  be  the  best  option.  Recall  that  in  particle  filtering 
we  start  with  N  i.i.d.  particles  distributed  according  to  the  initial  distribution.  In  the  resampling 
part  the  low  probability  particles  die  and  the  high  probability  particles  produce  many  particles 
identical  to  themselves.  Since  zt  does  not  change,  the  part  of  the  particles  associated  to  zt  tends 
to  cover  smaller  and  smaller  portions  of  the  state  space.  In  fact,  the  state  space  of  the  integer 
vectors  is  defined  by  the  particles  at  the  initial  time.  This  problem  can  be  overcome  by  modifying 
the  algorithm  mentioned  in  Section  (4).  In  the  new  algorithm,  Step  5  is  changed  in  such  a  way 
that  the  particles  are  the  addition  of  the  original  particles  found  by  Algorithm  4.1,  with  a  random 
vector.  The  modification  is  very  important  for  the  integer  values,  since  the  integers  do  not  have 
a  dynamics  that  is  driven  by  a  random  input.  In  [9],  a  similar  modification  has  been  used  for  the 
regular  nonlinear  filtering  setup.  It  seems  that  the  convergence  results  given  in  [9]  can  be  applied 
to  the  current  case  as  well. 

Based  on  the  modified  algorithm,  we  simulated  a  nonlinear  filtering  problem  similar  to  the 
problem  involved  in  the  GPS  system. 

In  a  two  dimensional  space,  three  transmitters  (imagine  three  pseudo  satellites)  are  mounted 
on  three  known  points  (2000,  100000),  (0,  100000),  and  (-2000,  100000).  The  moving  object  can 
measure  its  distance  from  these  transmitters.  For  each  pseudo  satellite,  two  types  of  measurement 
are  possible:  One  with  high  measurement  noise  and  the  other  with  low  measurement  noise.  For  the 
low  measurement  noise,  though,  there  is  an  integer  ambiguity.  The  dynamics  of  the  moving  object 
for  this  example  is  considered  to  be  in  discrete  time  and  linear  time  invariant.  The  dynamics  and 
observation  equation  is  given  as  follows: 

/  x i  \  /  1  At  0  0  \  /  x\  \  /  w\  \ 

Vl  _  0  1  0  0  Vl  VJ2 

X‘2  0  0  1  At  X2  103 

V  )n+1  V  0  0  0  1  /  \  ^/  „  V  )  n 

Vn  =||x-Si||+<  ,*  =  1,2,3 

ybn  =  llx-Si||  +  rii  +  vhn  ,  *  =  1,2,3, 

where  x  =  (xi,X2)T,  s *  is  the  position  of  pseudo  satellite  i  in  two  dimensional  space,  At  =  0.1 
unit  of  time,  n-i  is  the  integer  ambiguity  of  the  pseudo  satellite  i,  and  w  =  (wi,W2,W3,W4)T  and 

v  =  o\,  uf  j  *4,  uf ,  o^j  are  zero  mean  white  Gaussian  noise  process  with  covariance  matrices 

£w  =  drag  (1, 0.5, 1, 0.5)  and  Ev  =  diag  (5,0.2, 5,0.2, 5,0.2),  respectively.  In  the  simulation,  it  is 
assumed  that  the  initial  condition  for  the  position  is  distributed  in  a  square  of  size  200  x  200  units 
squared,  symmetric  with  respect  to  the  origin. 

In  brief,  the  simulation  can  be  separated  into  two  parts,  initialization  and  the  full  non-linear 
filtering.  In  the  initialization  part,  we  start  with  the  initial  probability  distribution  for  (xi,X2) 
and  from  a  series  of  observations,  we  find  an  estimate  for  the  probability  distribution  of  (01,02)- 
In  this  part,  we  do  not  use  the  dynamics  of  the  moving  object.  Using  our  estimate  for  the  prob¬ 
ability  distribution  of  (x\,  v\,  X2, 02)  we  find  the  distribution  for  the  integer  ambiguity.  After  this, 
the  initialization  is  over,  and  the  full  non-linear  filter  is  used.  There  are  some  minor  numerical 
considerations  that  we  would  like  to  point  out.  In  the  Bayes  step  of  the  algorithm,  the  numbers 
are  usually  very  small,  and  without  proper  scaling  the  original  algorithm  would  not  work.  In  the 
resampling  part,  one  can  use  the  law  of  large  numbers  and  regenerate  the  particles  based  on  their 
weight  without  generating  random  numbers  that  are  time  consuming.  The  result  of  the  simulations 
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The  Estimated  and  Actual  Integer  Ambiguity 


Figure  1:  Estimated  integer  ambiguity  versus  the  actual  integer  ambiguity  of  pseudo  satellite  (1). 
At  time  100  there  is  a  cycle  slip  of  strength  -20  for  the  measured  phase  of  the  carrier  from  pseudo 
satellite  (1). 


are  shown  in  Figures  1,  2,  3,  4,  5,  and  6.  To  display  the  estimated  integers,  we  simply  used  the 
mean  value,  which  is  not  necessarily  the  best  choice.  Of  course,  since  we  have  the  distribution, 
we  can  use  the  MAP  estimate  of  the  integers.  In  this  simulation  we  forced  one  of  the  integers  to 
have  a  jump.  Although  our  algorithm  is  not  designed  for  these  kinds  of  changes,  we  see  that  it  can 
estimate  the  new  integer  values.  In  future,  we  use  special  treatment  for  the  times  when  these  kinds 
of  jumps  happen.  As  we  can  see,  the  estimates  for  the  integers  are  reasonably  good.  The  reliability 
of  the  estimate  for  the  integers  depends  on  the  energy  of  the  noise. 


8  Future  Works 

The  simulations  results  show  that  our  method  is  capable  of  estimating  the  integer  ambiguity  and 
the  position.  There  are  certain  issues  that  need  further  investigation.  In  the  following,  we  itemize 
these  issues: 

•  What  are  the  proper  criteria  to  stop  the  integer  ambiguity  estimation  part  and  fix  the  integers? 

•  What  happens  when  a  cycle  slip  happens,  i.e.  one  or  more  of  the  integers  have  a  jump?  What 
change  detection  algorithm  is  proper  and  what  is  the  performance  of  this  algorithm?  How 
can  we  repair  the  integer  ambiguity  efficiently? 

•  What  happens  when  the  number  of  the  satellites  drops  from  the  critical  number? 

•  How  much  improvement  does  the  method  of  Section  5  for  integer  ambiguity  and  position 
estimation  have  over  particle  filtering? 
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Time 


Figure  2:  Estimated  integer  ambiguity  versus  the  actual  integer  ambiguity  of  pseudo  satellite  (2). 
At  time  100  there  is  a  cycle  slip  of  strength  -20  for  the  measured  phase  of  the  carrier  from  pseudo 
satellite  (1). 


Time 


Figure  3:  Estimated  integer  ambiguity  versus  the  actual  integer  ambiguity  of  pseudo  satellite  (3). 
At  time  100  there  is  a  cycle  slip  of  strength  -20  for  the  measured  phase  of  the  carrier  from  pseudo 
satellite  (1). 


24 


Th  estimated  x  component  versus  the  actual  x  component 


Figure  4:  Estimated  x\  component  versus  the  actual  x\  component  of  the  position  of  the  car.  At 
time  100  there  is  a  cycle  slip  of  strength  -20  for  the  measured  phase  of  the  carrier  from  pseudo 
satellite  (1). 


The  estimated  y  component  versus  the  actual  y  component 


Figure  5:  Estimated  X2  component  versus  the  actual  X2  component  of  the  position  of  the  car.  At 
time  100  there  is  a  cycle  slip  of  strength  -20  for  the  measured  phase  of  the  carrier  from  pseudo 
satellite  (1). 
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The  actual  trajectory  versus  the  estimated  trajectory 


Figure  6:  Estimated  trajectory  versus  the  actual  trajectory  of  the  car.  At  time  100  there  is  a  cycle 
slip  of  strength  -20  for  the  measured  phase  of  the  carrier  from  pseudo  satellite  (1). 


These  questions  are  to  be  answered  in  the  future  work.  In  addition  to  these,  we  shall  be  more 
specific  in  our  simulations,  and  use  real  GPS  data  for  our  results. 
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